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The spreading of polymer droplets is studied using molecular dynamics simulations. To study the 
dynamics of both the precursor foot and the bulk droplet, large drops of 200,000 monomers are 
simulated using a bead-spring model for polymers of chain length 10, 20, and 40 monomers per chain. 
We compare spreading on flat and atomistic surfaces, chain length effects, and different applications 
of the Langevin and dissipative particle dynamics thermostats. We find diffusive behavior for the 
precursor foot and good agreement with the molecular kinetic model of droplet spreading using 
both flat and atomistic surfaces. Despite the large system size and long simulation time relative to 
previous simulations, we find no evidence of hydrodynamic behavior in the spreading droplet. 
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I. INTRODUCTION 

The spreading of liquid droplets on a surface is an im- 
portant issue for several industries including adhesion, 
lubrication, coating, and printing. Emerging nanotech- 
nology in areas such as lithography and microfluidics has 
made the issue of droplet spreading on small length scales 
even more relevant. Experiments on droplet spreading 
have revealed several phenomena involved in the spread- 
ing process, some of which occur on the atomic level and 
others that become relevant at meso scop ic length scales 
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include the spreading of a precursor foot ahead of the 
dr oplet El , terraced spreading of mono- molecular layers 
[il llrllrn f. and viscous losses due to rolling motion [Ulisj. 

Several models have been proposed to describe the 
spontaneous spreading of liquid droplets on a surface. 
These models can be classified as molecular kinetic mod- 
els, continuum hydrodynamic models, or combined mod- 
els. The molecular kinetic theory of Eyring [19( has been 
applied to the kinetics of wetting by Blake and Haynes 
[2fl E3 as well as Cherry and Holmes H^. This the- 
ory treats the surface adsorption of liquid molecules as 
the dominant factor in the spreading of a droplet. The 
hydrodynamic theory 0, 0, [2^, 0] focuses on the en- 
ergy dissipation due to viscous flow in the droplet. It 
has been claimed that hydrodynamic dissipation is dom- 
inant for small contact angles and non-hydrodynamic dis- 
sipation is dominant for relatively large contact angles 
|25|. Since both mechanisms are present in spreading 
droplets, several grou ps have proposed combined theo- 
ries 0, H, l2rj l27l |2S| . Experimental results for the 
spreading of poly(dimethylsiloxane) (PDMS) drops on 
bare silicon wafers have shown good agreement with one 
combined model jToj | . 

The study of droplet spreading using molecular dy- 
namics simulation has been hindered due to computa- 
tional limitations restricting simulations to small droplet 
sizes and short times. Molecular dynamics simulations 
were first used to study the spreading of monomer and 
dimer liquids J2jj EiJ, Ell However, the spreading 

of monomer and dimer droplets are clearly influenced by 



the volatility of the small molecules, allowing them to 
vaporize and condense independent of the dynamics of 
the droplet. To separate the spreading from the vapor- 
ization and condensation, subsequent simulations used 
short bead-spring chain molecules since they have a very 
low vapor pressure. In most cases, the simulations re- 
produced the experimentally observed R ~ t 1 / 2 scaling 
of the contact radius of the precursor foot on both atom- 
istic [lfl |H H3 and flat |M E3 surfaces, though 
logarithmic scaling has also been observed j3f|. It is be- 
lieved that this difference is due to the corrugation of the 
substrate, producing i 1 / 2 scaling for a sufficiently small 
lattice dimension and a logarithmic scaling for larg e, i.e. 
strong corrugation 38]. Milchev and Binder 39] have 
studied wetting using Monte Carlo simulations on a flat 
substrate which suggest Tanner's spreading law for the 
growth dynamics of the droplet holds on the nanoscopic 
scale. Other comparisons to theoretical models have 
strongly suppor ted the molecular kinetic theory of wet- 
ting O El |43|, Q, probably due to the relatively 
small droplet sizes and short simulation times employed. 

In this paper, we present results from extensive molec- 
ular dynamics simulations of coarse-grained models of 
polymer droplets wetting a surface. Although most re- 
cent simulations of droplet spreadin g us e droplets con- 
taining 20, 000 to 32, 000 monomers [H |M H El 113 
E2, we consider drops composed of 100,000 to 200,000 
monomers to simultaneously study the precursor foot and 
bulk regions for long times. We compare simulations 
performed using both a flat surface and an atomistic 
substrate to determine if the computationally expensive 
atomistic substrate is required to obtain correct spread- 
ing dynamics. We also evaluate different implementa- 
tions of the Langevin and dissipative particle dynamics 
(DPD) thermostats for efficiency and realism in preserv- 
ing hydrodynamic effects. Also, the difference in using a 
spherical droplet as the starting configuration as opposed 
to a hemispherical droplet is discussed. We find that the 
method which captures all of the physics of the spreading 
drop in the most computationally efficient manner is to 
simulate large drops on flat substrates with a coupling to 
the thermostat which falls off exponentially with distance 
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from the substrate [45j . For atomistic substrates, we find 
coupling only the substrate monomers to the Langevin 
thermostat significantly more efficient than coupling the 
DPD thermostat to all monomers. 

The paper is organized as follows. Section ITT1 describes 
the details of the molecular dynamics simulations and 
the application of the thermostats. Section HTD presents 
the results for the time dependence of the contact ra- 
dius. The contact angle data is fit to models of droplet 
spreading in Section llVl and conclusions are presented in 
Section E| 



II. SIMULATION DETAILS 
A. System 

We perform molecular dynamics (MD) simulations us- 
ing a coarse-grained model for the polymer chains in 
which the polymer is represented by spherical beads of 
mass m attached by springs. We use a cutoff Lennard- 
Jones (LJ) potential to describe the interaction between 
all monomers. The LJ potential is given by 
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where e a p and a a p are the LJ units of energy and length 
and the cutoff is set to r c = 2.5 a a p. We denote the 
polymer monomers as type 1 and substrate monomers as 
type 2. The monomer-monomer interaction, e±i = e, is 
used as the reference and all monomers have the same 
diameter a a p = a. For bonded monomers, we apply an 
additional potential where each bond is described by the 
finite extensible nonlinear elastic (FENE) potential |4(| . 



the wetting droplets, we study two substrates, contain- 
ing either 49 200 or 99 960 monomers per layer. The di- 
mensions of the substrates are 23f.2 a x 23f .0 ct and 
330.8cr x 33f.4cr, respectively. We refer to these as the 
small, medium and large substrates. The large substrates 
are necessary because the finite size of the atomistic sub- 
strates require the use of periodic boundary conditions 
at their edges whereas the flat surface can extend in- 
definitely in the x and y directions. For the atomistic 
substrate, during the course of the simulation, the pre- 
cursor foot reaches the edge of the substrate and interacts 
with the periodic image of the droplet. Although this can 
be related to the spreading of an array of nanodroplets, 
such as in micro-contact printing, we do not include any 
data for the precursor foot once it reaches the periodic 
image. The droplets consist of ~ 100,000 monomers for 
non-wetting droplets and ~ 200, 000 monomers for wet- 
ting droplets. All simulations are run at a temperature 
of T= 1.0e/k B - 

For the flat surface, the interaction between the 
monomers in the droplet and the surface is modeled by 
an integrated LJ potential, 
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with z c — 2.2cr. 

The equations of motion are integrated using a 
velocity- Verlet algorithm. We use a time step of At = 

1/2 

0.009 r where t — a (^) • The simulations are per- 
formed using the lammps code on 36 to 100 Dec 
Alpha processors of Sandia's C Plant cluster. Simulating 
one million steps for a wetting drop of 200, 000 monomers 
on the medium atomistic substrate takes between 90 and 
250 hours on 64 processors, depending on the thermostat. 
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with k = 30 e and Rq — 1.5 ct. 

Droplets consisting of chains of length N = 10, 20, or 
40 monomers per chain are created by first equilibrat- 
ing a melt of the polymer and then removing molecules 
whose centers are outside of a hemisphere of a given ra- 
dius, 38 cr for non- wetting droplets and 48 ct for wetting 
droplets. The droplet is then placed on either an atom- 
istic substrate or a flat substrate. 

The atomistic substrate is composed of LJ particles 
forming four layers of the (111) surface of an fee lattice 
where the bottom layer is frozen and the top three layers 
maintain their structure through a strong LJ interaction, 
£22 = 5e. The masses of the substrate monomers are set 
to m,2 = 2mi = 2m. For non-wetting droplets, each 
layer of the substrate contains 12 000 monomers and the 
dimensions of the substrate are 110. 0(7 x 115. 4ct. For 



B. Thermostats 

The choice of thermostat employed can greatly affect 
the droplet spreading dynamics, so we compare simu- 
lations that use the Langevin 48| and DPD [ZiJ H3 
thermostats. The purpose is to find an approach that 
is both computationally efficient and provides a realistic 
representation of the transfer of energy in the spreading 
droplet. 

The Langevin thermostat simulates a heat bath by 
adding Gaussian white noise and friction terms to the 
equation of motion, 



mrfi = -AUi - mcijjii + Wi(t), 



(4) 



where 7^ is the friction parameter for the Langevin ther- 
mostat, —AUi is ft 16 force acting on monomer i due to 
the potentials defined above, and Wj(t) is a Gaussian 
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white noise term such that 



(Wi(i) • Wj(t')) = &k B Tm tlL 5 t] 5 (t - t') . 
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The Langevin thermostat can either be coupled to all 
monomers in the system or just to those in the sub- 
strate. The advantage of the latter is that the long-range 
hydrodynamic interactions are preserved in the droplet, 
whereas coupling all monomers to the Langevin thermo- 
stat screens the hydrodynamic interactions. Both ap- 
proaches are applied in the simulations to test the various 
models for droplet spreading discussed below in Sec. IV. 
The damping constant is chosen to be jl — 0.1 r _1 in 
most cases, which is much smaller than that arising from 
collisions between monomers. 

Our next approach is to apply the thermostat from the 
DPD simulation method. The DPD technique includes 
a dissipative force term in the equations of motion along 
with random forces. The equation of motion for the DPD 
thermostat is 



m l v l = (-&Uij + + F 
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In Eq. |SJ F^ and F^ are the dissipative and random 
terms given by 



FPj = -ma D PDW 2 (r i3 ) (f y - • (r* - tj)) fy (7) 



m i °'DPD'w(rij)Ci : 



(8) 



where jdpd is the DPD friction parameter, <Jpp D — 
2kBT"fDPD, Qj is a Gaussian noise term with 
(Qj{tXki{t')) = (SikSji+SuSj^Sit-t), r y - - rt-rj, 
r ij/ r ij- The weight function w(rij) 
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We take r' c = r c = 2.5 a. The advantage of this ther- 
mostatting technique is that the momentum is conserved 
locally and long-range hydrodynamic interactions are 
preserved even in the case where all monomers are cou- 
pled to the thermostat. All simulations with the DPD 
thermostat use jdpd = 0.1 r _1 , so the dissipation from 
the thermostat is much less than from monomer collisions 
as seen in Sec. IIVI Simulations that use DPD couple the 
thermostat to all atoms in the system. 

In the case of the flat substrate, we study several 
methods to thermostat the system. In the first, we 
simply couple the Langevin or DPD thermostat to all 
monomers. However this is somewhat unphysical since 
monomers near the substrate are expected to have a 



stronger damping than those in the bulk of the droplet. 
In the case of the Langevin thermostat, this coupling of 
all monomers also means that the hydrodynamic inter- 
actions are screened. In addition, chains which separate 
from the droplet move across the substrate very rapidly, 
particularly for the DPD thermostat. For this reason, we 
did not further pursue the DPD thermostat on the flat 
substrate. To overcome these difficulties, we follow the 
approach of Braun and Peynard |45| and add an external 
Langevin coupling with a damping rate that decreases 
exponentially away from the substrate. We choose the 
form 



7l(z) = 7£exp(cr - z) 



(10) 



where 7£ is the surface Langevin coupling and z is 
the distance from the substrate. We choose values of 
7£ = 1.0, 3.0, and 10.0 r -1 . There is no obvious a priori 
way to define the appropriate value of 7£. However, one 
way is to choose 7£ so that the diffusion constant of the 
precursor foot is comparable for the flat and atomic sub- 
strates for comparable departures from the wetting/non- 
wetting transition (see Fig. 0] below). 



III. SIMULATION RESULTS 

A droplet containing about 200, 000 monomers for a 
wetting droplet is large enough to allow us to simulta- 
neously study the bulk and precursor foot regions. This 
can be seen in the profile views for chain length N = 10 
in Figs. Hand [3 which show the foot extending beyond 
the bulk region for wetting droplets on an atomistic sub- 
strate and a flat surface, respectively. Note that Fig. Q] 
shows the thickness of the foot increasing after it reaches 
the periodic image. The same behavior is seen when pe- 
riodic boundaries are applied to the flat surface, so this 
is not an effect of the corrugation of the substrate. 

To characterize the spreading dynamics of these 
droplets, we extract the instantaneous contact radius and 
contact angle every 10, 000 to 40, 000 At. The contact 
radius is calculated by defining a two-dimensional radial 
distribution function, g(r) — p{r)/ p 1 based on every par- 
ticle within 1.5 a of the surface. The local density at a 
distance r from the center of mass of the droplet is 
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where N(r) is the number of particles at a distance be- 
tween r and r + Ar from the center of mass and p is 
the integral of p(r) over the entire surface. The contact 
radius is defined as the distance r at which g(r) = 0.98. 
This approach provides a robust measure of the radius 
at any point during the spreading simulation. The same 
calculation is used to obtain the droplet radius for ten 
slices of the droplet at incremental heights every 1.5 a 
from the surface. A line is fit to the resulting points and 
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FIG. 1: Profile of the N = 10 polymer droplet spreading 
on the atomistic substrate at three different times using the 
Langevin thermostat applied only to the substrate monomers 
with 7l = 0.1 r -1 and E12 = 1.5 e. 
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FIG. 2: Profile of the N = 10 polymer droplet spreading 
on the flat surface at three different times using the surface 
Langevin thermostat. 7£ = 10.0 t~ , e w = 2.0 e. 



the instantaneous contact angle is determined from the 
slope of the line. For simulations that exhibit a precursor 
foot, the particles within 4.5 a of the surface are ignored 
in the contact angle calculation. 

The non-wetting droplets reach their equilibrium con- 
figurations fairly rapidly, as shown by the contact angle 
data in Fig. [3J The equilibrium contact angle measured 
as a function of polymer-surface interaction strength is 
shown in Fig. 0] From this figure, it is clear that 
the transition from non-wetting to wetting occurs near 
e\ 2 — 1-05 e for droplets on a substrate and ~ 1.75 e 
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FIG. 3: Contact angle of non-wetting droplets of N = 10 poly- 
mers on an atomistic substrate starting from a hemispherical 
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FIG. 4: Equilibrium contact angle as a function of polymer- 
surface interaction strength showing the transition from non- 
wetting to wetting for N = 10 polymer droplets on an atom- 
istic substrate (circles) and a flat surface (squares) . 



for droplets on a flat surface. For most of the wetting 
simulations, we use £12 = 1.5 £ for the atomistic sub- 
strate and e w = 2.0 e for the flat substrate, both well 
within the wetting regime. 

The time dependence of the contact radius of the pre- 
cursor foot and bulk region is shown in Fig. [5]for wetting 
droplets on an atomistic substrate for three chain lengths. 
The i 1 / 2 behavior is evident for the precursor foot at 
all chain sizes, while the kinetics of the main droplet is 
clearly significantly slower. The ./V = 10 data shown in 
Fig. is taken from simulations on both the large and 
medium substrates whereas the TV = 20 and N = 40 sim- 
ulations are on the medium substrate. The contact radius 
of the bulk droplet increases steadily for all three chain 
lengths on the medium substrate. However, the run on 
the large substrate shows a slowing down and eventual 
contraction of the bulk contact radius as the foot con- 




FIG. 5: Time dependence of the contact radius of the pre- 
cursor foot and bulk droplet for wetting droplets on an atom- 
istic substrate at three different chain lengths starting from a 
hemisphere with a contact angle of approximately 90°. The 
Langevin thermostat, jl = O.lr -1 , is applied only to the sub- 
strate monomers and £12 = 1.5 e. Results for N = 10 are for 
both the medium and large atomistic substrates, while those 
for N = 20 and 40 are for the medium atomistic substrate. 
The inset shows the contact radius for N = 10 starting with 
a spherical droplet (solid line) compared to a hemispherical 
droplet (dotted). Results for the hemisphere in the inset have 
been shifted downward to easily compare the late time behav- 
ior. 

tinues outward, depleting the supply of material in the 
bulk faster than the drop can transfer material down- 
ward. This suggests that for our largest substrate, the 
drop size must be even larger to be able to study both the 
precursor foot and bulk droplet in the same simulation. 

The inset in Fig. \5\ shows the spreading of a spherical 
N = 10 droplet compared to an initial hemisphere. The 
sphere is placed just above the substrate with zero initial 
velocity to avoid any effect due to impact velocity. The 
difficulty in measuring the spreading rate for this case 
is evident as it takes roughly 1200 r for the sphere to 
adopt a hemispherical shape, 1600 r for the spreading 
rate of the foot to match that of the hemisphere, and 
5000 r for the spreading rate of the bulk to match that of 
the hemisphere. (The hemisphere data for the foot and 
bulk regions are shifted downward to easily compare the 
spreading rates.) 

Voue et al. 0, El found both experimentally for 
PDMS droplets and in numerical computer simulations 
that the diffusion constant of the precursor foot varies 
non-monotonically with increasing coupling to the sub- 
strate. At first, increasing the coupling to the substrate 
increases the driving force and the fluid spreads on the 
substrate more rapidly. However, further increases in 
the strength of the fluid substrate coupling, while in- 
creasing the driving force, also increase the friction of 
the fluid monomers with the substrate, resulting in a de- 
crease in the diffusion constant. From the time depen- 
dence of R(t), the diffusion constant Df for the foot can 



t 1 r 




1 1 1 1 1 1 1 

50 100 150 200 



FIG. 6: Effect of thermostat on contact radius of precursor 
foot and bulk region for wetting droplets of N = 10 polymers 
on an atomistic substrate for £12 = 1.5 £. The thermostats 
applied are DPD (solid line) , Langevin on all monomers (dot- 
ted) and Langevin on only substrate monomers (dashed). 

JDPD =71, = 0.1 T -1 . 

be determined from 

((JZ(t) - i?(Q)) 2 ) = AD jt (12) 

The resulting diffusion constants are Dt = 0.34 <t 2 /t, 
0.30 ct 2 /t, and 0.23 cr 2 /r for N = 10, 20, and 40 respec- 
tively for £12 = 1.5 e, indicating a very weak dependence 
on chain length, at least for these unentangled chains. In- 
creasing £12 to 2.0 e, we find Df — 0.16 <j 2 /t for N = 10, 
thus the droplets are in the high friction regime for these 
values of fluid substrate coupling. 

Figure shows the time dependence of the contact ra- 
dius for wetting droplets on an atomistic substrate using 
different thermostatting techniques. These results show 
that there is essentially no difference in the spreading rate 
between the DPD thermostat applied to all monomers 
and the Langevin thermostat applied only to the sub- 
strate. We can see that applying the Langevin ther- 
mostat to all monomers slightly decreases the spreading 
rate as the viscous heating is removed from the system, 
though the resulting loss of hydrodynamic flow, at least 
for the droplet size studied here, has no significant im- 
pact. 

For wetting droplets on a flat surface, the thermostat 
dependence of the contact radius is shown in Fig. [7| 
Here, the Langevin thermostat is applied either to all 
monomers (curves labeled with jl) or with the surface 
Langevin coupling (curves labeled with 7J). The value of 
7l clearly has a strong influence on the spreading rate. 
7£ = 3.0 r^ 1 gives a diffusion constant comparable to 
the atomistic substrate with 7^ = 0.1 t _1 . The chain 
length dependence of the contact radius is shown in Fig. 
[S] Again, the t 1 / 2 behavior is evident in the foot region 
but not the bulk region. The chain length dependence 
on the flat surface is similar to the atomistic substrate, 
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FIG. 7: Effect of thermostat on contact radius of (a) precur- 
sor foot and (b) bulk region for wetting droplets of N = 10 
polymers on a flat surface with s w = 2.0 e. 



showing a moderate decrease in spreading rate for larger 
polymers. 



IV. MODELS OF DROPLET SPREADING 
DYNAMICS 

A. Overview of models 

The dynamics of droplet spreading are controlled by 
the driving force (the difference in surface tension 7 at 
each interface) and by the energy dissipation. The to- 
tal energy dissipation can be represented by a sum of 

three different components, T (t w +t f +ti^ 3]. The 

first term, TT, W , represents energy dissipation due to the 
hydrodynamic flow in the bulk of the droplet as more 
material is transferred to the surface. TE f relates to the 
viscous dissipation in the precursor foot present in cases 
of complete wetting. The third term, TE/, refers to the 
dissipation in the vicinity of the contact line due to the 
adsorption and desorption of liquid molecules to the solid 
surface. Here, we compare models that incorporate one 
or more of these dissipation mechanisms to our simula- 
tion results. 
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FIG. 8: Chain length dependence of the contact radius of the 
precursor foot and bulk droplet for wetting droplets on a flat 
surface with e w = 2.0 e. The surface Langevin thermostat is 
applied with = 10.0 r -1 . 



The molecular kinetic theory of liquids developed by 
Eyring and coworkers 0] has been applied to droplet 
spreading by Blake and Haynes [2]]]. It focuses on the 
adsorption of liquid molecules to the surface as the dom- 
inant factor in energy dissipation. In this theory, the 
liquid molecules jump between surface sites separated by 
a distance A with a frequency K . The velocity of the 
contact line is related to the contact angle 9 by 



dR 

— = 2i4TAsinh 
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where 7 is the surface tension of the liquid/vapor inter- 
face, An is the density of sites on the solid surface, and 
9q is the equilibrium contact angle. For sufficiently low 
velocities, the equation can be written in its linearized 
form, 
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Assuming the droplet maintains constant volume and the 
shape of a spherical cap, the velocity of the contact line 
can be expressed in terms of the time dependence of the 
contact angle purely from geometric arguments giving 
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Combining Eqs. 1141 and 1151 gives an expression for the 
time dependence of the contact angle, 

f = - (^7) V3 O (9) f Q (cos 9 - cos9) (16) 



where 
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and Co is the friction coefficient defined as Co 
which has units of viscosity. 

The hydrodynamic model [24j describes the flow pat- 
tern that forms in the bulk of the droplet as material is 
transferred to the advancing contact line. This model 
can be obtained by solving the equations of motion and 
continuity for the droplet described as a cylindrical disk 
[5ll ] instead of a spherical cap. Neglecting the flow per- 
pendicular to the surface and balancing the radial shear 
stress at the top of the cylinder with the effective radial 
surface tension, the velocity of the contact line is written 
as 
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where V is the droplet volume, r\ is the viscosity of the 
liquid, and (3=1 — cos#o- Equation 1181 is in agree- 
ment with Tanner's spreading law for completely wet- 
ting systems (#0 = 0) and for non-wetting systems with 
small equilibrium contact angles, giving R ~ t 1 ' 10 at long 
times. Instead of directly combining Eqs. IT31and[THl we 
apply the approach of de Ruijter et al. [g, |l(j in order 
to make a direct comparison with the combined model 
presented below. Using the same cylindrical disk model, 
they neglect the flow perpendicular to the surface and 
specify that the velocity at the upper edge of the cylin- 
der is the actual droplet spreading rate, dR/dt. With 
this approach, they find that the hydrodynamic dissipa- 
tion term can be written as 
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where <f> (#)is a geometric factor defined as 
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and a is an adjustable parameter that represents the ra- 
dius of the core region of the droplet, where the radial 
flow is negligible. For the hydrodynamic model, they 
obtain 



dt 



= -( 
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Both types of dissipation are present in the spreading 
droplet. The hydrodynamic mechanism is expected to 
dominate at low velocities and small contact angles while 
the kinetic mechanism is expected to dominate at high 
velocities and large contact angles [25|. We include in 
our comparison a model developed by de Ruijter et al. 
0, ITol ] containing both kinetic and hydrodynamic terms. 
In this model, the velocity of the contact line is written 
as 



dR 7 [cos 9 — cos 6] 
~dt ~ Co + 6770 (0) In [R/a] 
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TABLE I: Bulk properties of bead-spring chains obtained 
from MD simulation for T = e/fcg, P — 0. 
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0.92±0.02 


17.4±0.7 


3.04±0.03 


16.4 


40 


0.8856 


0.95±0.02 


41.7±1.4 


1.23±0.01 


20.4 



Combining this with Eq. 1151 gives 



d9 r tt \ 



dt 



\3VJ 



1 ^ 3 r,tn\ J {COS 9 — COS 8) 

Co + 6770(0) In [R/a] 



0,(0) 



(23) 



B. Analysis of Models 

Fitting simulation data to the models described above 
requires both the liquid/ vapor surface tension and the 
bulk viscosity of the polymer. The surface tension, 7, 
is obtained by first constructing a slab of the polymer 
melt containing 10, 000 chains of N = 10, 5000 chains of 
N = 20, or 5000 chains of N — 40 centered in the simu- 
lation box such that there are two surfaces perpendicular 
to the z direction. The simulations are run at tempera- 
ture T = 1.0 and pressure P ~ without tail corrections. 
We leave out tail corrections to the pressure in order to 
match the system of the spreading droplet. The simu- 
lations are run until the two liquid/ vapor interfaces are 
equilibrated, as determined by the density profiles across 
the interfaces. From the equilibrium values of the pres- 
sure, parallel and perpendicular to the interfaces, 7 can 
easily be determined from |52j 



7 



L, 



[p± ( z ) ~P\\ (z)] dz. 



(24) 



The values for the surface tension are summarized in Ta- 
ble||J These values can be compared to 7 = 0.08e/cr 2 for 
a system of monomers |53|. 

The viscosity is computed from the equilibrium fluc- 
tuations of the off-diagonal components of the pressure 
tensor |54| . The pressure tensors are recorded from sim- 
ulations of systems containing melts of 500 chains of the 
N = 10 polymer, 250 chains of the N — 20 polymer, 
and 500 chains of the N = 40 polymer at T = 1.0 
with the bulk pressure P ~ without tail corrections. 
These simulations are run at a timestep At = 0.006 for 
up to 25, 000 r. The autocorrelation function of each 
off-diagonal component of the stress tensor is calculated 
using the Numerical Recipes routine CORREL [55]. The 
autocorrelation functions are averaged to improve statis- 
tical uncertainty. From this, the viscosity can be calcu- 
lated using ,54] 



V 



V = 



dt (<7 a p{t)(T a p{0)) 



ksT Jo 

The results for r\ are summarized in Table [I] 



(25) 
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+ y =10.0 x N=40 



7, =10.0 x N=20 
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FIG. 9: Fits to contact angle data (symbols) of (a) kinetic, (b) 
hydro dynamic and (c) combined models for wetting droplets 
on a flat surface with e w = 2.0 e. The chain length is N — 10 
unless otherwise specified. The Langevin thermostat is ap- 
plied to all monomers (jl) or just monomers near the surface 
(7!,). The data sets are shifted by 10° increments (except for 
Jl = 1.0 r" 1 ) for clarity. 



Estimates of the friction coefficients, Cr, obtained from 
the melt simulations, are included in Table UJ The diffu- 
sion constant D is determined from the mean square dis- 
placement of the middle monomers of each chain and us- 
ing the Rouse model one can extract £r from D = 1 ' 



mN(a 



With the above values for the surface tension and vis- 
cosity, the simulation data is fit to each of the models de- 
scribed above. The fit is performed by taking initial guess 
values for the independent parameters and integrating 
the expression for d9/dt defined in one of the equations 
Hoi [?T1 or The integration uses the fourth-order 

Runge-Kutta method to generate a set of data, 9 ca i c (t). 
The parameters are varied using the downhill simplex 
method [55| until the difference between the model and 
simulation data, \0 ca i c {t) — 6(t)\ /0(t), is minimized. 
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FIG. 10: Fits to contact angle data (symbols) of (a) ki- 
netic, (b) hydrodynamic and (c) combined models for wetting 
droplets on the medium atomistic substrate with £12 = 1.5 e. 
The Langevin thermostat is applied to all monomers or to 
just substrate monomers. 71, = 0.1 t -1 for all cases except 
DPD where ~/dpd = 0.1 r" 1 . The data sets are shifted by 
10° increments (except for N — 10 DPD) for clarity. 



The kinetic, hydrodynamic and combined models are 
fit to the contact angles of droplets spreading on a flat 
surface in Fig. |5J The Langevin thermostat is applied 
either to all monomers or only to those near the surface. 
We find that both the kinetic and combined models fit the 
data well despite the fact that they predict that the fric- 
tion coefficient, £oj is larger in the combined model than 
in the kinetic model. The hydrodynamic model produces 
a very poor fit to each data set as shown in Fig. [SJx The 
best fit parameters for these models applied to data for 
wetting droplets on a flat surface are shown in Table ITU 
The error reported for each model is calculated as 



X 



1 

A7 



M 

£ 

i=l 



,Ut)-o{t)Y 

8(t) 



(26) 



where Af is the number of data points in each set of data. 

Figure HUl shows the kinetic, hydrodynamic, and com- 
bined model fits to the contact angle data for wetting 
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TABLE II: Model parameters and error estimates resulting from fits to contact angle data from simulations of wetting droplets 
on a fiat surface. Values for 7^ and 7£ are listed in the first column. e w = 2.0 e. 







Kinetic 


Hydrodynamic 


Combined 








Thermostat 


N 


Co (£) 


a(er) 


Co (7^) 


a(cr) 


Xlin 


Xhydro 


Xcomb 


7l = 1.0 r- 1 


10 


9.55 


44.40 


35.41 


71.23 


0.0022 


0.047 


0.0039 


7£ = 3.0r- 1 


10 


25.97 


42.29 


57.27 


71.55 


0.00028 


0.025 


0.0011 


7£ = 10.0 r- 1 


10 


56.30 


38.14 


89.98 


83.80 


0.00024 


0.018 


0.00028 


7£ = 10.0 r- L 


20 


81.37 


38.83 


137.99 


84.77 


0.00015 


0.015 


0.00023 


Yl = 10.0 r^ 1 


40 


101.29 


38.63 


200.45 


86.49 


0.00022 


0.024 


0.00036 



droplets on the medium substrate. Again, the hydrody- 
namic model gives a significantly worse fit to the data. 
The best fit parameters for these models for wetting 
droplets on an atomistic substrate are shown in Table 
IIIII The kinetic and combined models give friction coef- 
ficients that are generally larger than the bulk viscosity 
for the range of coupling parameters used here. We find 
that, in contrast with previous work by De Ruijter et al. 
0, El , the combined model predicts a larger friction co- 
efficient than the kinetic model. Also, the hydrodynamic 
and combined models give a value of a that is on the 
order of the radius of the droplet, indicating that hydro- 
dynamic flow is not a dominant feature of the spreading 
of these droplets, at least for the time scales accessible to 
simulation. 



sonable values for the friction coefficients, which we veri- 
fied through separate polymer melt simulations. Using a 
combined model that adds a hydrodynamic energy dissi- 
pation mechanism slightly improves the fit, but resulted 
in less accurate estimates of the friction coefficients. The 
fact that we do not observe evidence of hydrodynamic 
flow behavior may be due to the small droplet sizes ac- 
cessible to molecular dynamics simulation. Evidence for 
hydrodynamic effects on spreading has been observed 
experimentally for macroscopic drops 0, 0, H3- The 
length scale where hydrodynamic effects become impor- 
tant remains an open question. 

Future work will include studying the spreading be- 
havior of binary droplets and developing more realistic 
surface interactions. 



CONCLUSIONS 



VI. ACKNOWLEDGEMENTS 



In this study, we perform molecular dynamics simula- 
tions of polymer droplets that are roughly an order of 
magnitude greater in size than those previously stud- 
ied. We find this to be necessary to adequately model 
the behavior of the precursor foot and the bulk material 
simultaneously. Starting from a hemispherical droplet, 
we find that the precursor foot forms immediately and 
spreads diffusively for each system where the surface in- 
teraction strength is above the wetting/non- wetting tran- 
sition. The bulk region of the droplet spreads at a sig- 
nificantly slower rate, but the data is too imprecise to 
distinguish between, for example, a t 1 / 7 or a i 1 / 10 scaling. 

We perform spreading simulations on both an atom- 
istically realistic substrate and a perfectly flat surface. 
The simulations using a flat surface exhibit the same be- 
havior as the realistic substrate and greatly improve the 
computational efficiency since the number of monomers 
on the realistic substrate is typically several times greater 
than the number of monomers in the droplet. However, 
to do so, it is critical to apply a thermostat that cou- 
ples only to monomers near the surface. On an atomistic 
substrate, the most efficient method is to couple only the 
substrate particles to the thermostat. This is computa- 
tionally faster than coupling all monomers to the DPD 
thermostat and leads to the same results. 

Several droplet spreading models have been developed 
to fit contact angle data. A simple kinetic mechanism 
for energy dissipation fits the data well and provides rea- 
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dia is a multiprogram laboratory operated by Sandia Cor- 
poration, a Lockheed Martin Company, for the United 
States Department of Energy's National Nuclear Se- 
curity Administration under Contract No. DE-AC04- 
94AL85000. 
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TABLE III: 

Model parameters and error estimates resulting from fits to contact angle data from simulations of wetting droplets on an 



atomistic substrate. £12 = 1.5 e, 7dpd = 7l = 0.1 r 







Kinetic 


Hydrodynamic 


Combined 








Thermostat 


N 


Co (£) 


a(cr) 


Co (£) 


a(o-) 


xl in 


/(.hydro 


/Ccomb 


DPD 


10 


36.7 


42.1 


53.4 


65.5 


0.001 


0.015 


0.001 


Lang on All 


10 


50.8 


38.1 


81.8 


70.0 


0.001 


0.015 


0.001 


Lang on Sub 


10 


38.0 


41.8 


64.9 


69.6 


0.001 


0.020 


0.001 


Lang on Sub 


20 


54.4 


43.2 


91.9 


68.3 


0.001 


0.020 


0.001 


Lang on Sub 


40 


65.5 


42.9 


126 


64.0 


0.002 


0.019 


0.002 
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